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ABSTRACT 

Context. High-contrast imaging for the detection and characterization of exoplanets relies on the instrument's capability 
to block out the light of the host star. Some current post-processing methods for calibrating out the residual speckles 
use information redundancy offered by multispectral imaging but do not use any prior information on the origin of 
these speckles. 

Aims. We investigate whether additional information on the system and image formation process can be used to more 
finely exploit the multispectral information. 

Methods. We developed an inversion method in a Bayesian framework that is based on an analytical imaging model to 
estimate both the speckles and the object map. The model links the instrumental aberrations to the speckle pattern in 
the image focal plane, distinguishing between aberrations upstream and downstream of the coronagraph. 
Results. We propose and validate several numerical techniques to handle the difficult minimization problems of phase 
retrieval and achieve a contrast of 10 6 at 0.2 arcsec from simulated images, in the presence of photon noise. 
Conclusions. This opens up the the possibility of tests on real data where the ultimate performance may override the 
current techniques if the instrument has good and stable coronagraphic imaging quality. This paves the way for new 
astrophysical exploitations or even new designs for future instruments. 

Key words. High angular resolution - Image processing - Detection 



1. Introduction 

Ground-based instruments have now demonstrated the 
capability of detec t ing planetary mass companions 
(IChauvin et all l2004t lLagrange etall 120101: iMarois et all 
120011 " around bright host stars. By combining adaptive op- 
tics (AO) system and coronagraphs, some first direct de- 
tections from the ground have been possible in favorable 
cases, at large separations and in young systems when low- 
mass companions are still warm (> 1000 K) and therefore 
not too faint. There is a very strong astrophysical case to 
improve the high-contrast detection capability (10 5 for a 
young giant planet to 10 10 for an earth-like planet in the 
near infrared) very close to stars (< 0.1" to 1"). 

Several instruments will be capable of perform- 
ing multispectral imaging and will allow character- 
izing the planets by measuring their spectra . Thi s 
is the case of GPI ( Gemini) dGraham et all [2007h . 
Palm 3000 (Palomarl dHinklev et all 1201 lib SCExAO 
fSubaru) (iMartinache fc Guvonl I2009D . SPHERE 
(VLT) (jBeuzit et all I2008D . and several others tha t 
will follow, such as EPICS (E-ELT) (iKasper et al.ll2008lh 
By combining extreme adaptive optics (Ex-AO) and 
more accurate coronagraphs than before, the level of star 
light cancellation is highly improved, leading to a better 
signal-to-noise ratio. Even so, the residual host star light 
is affected by the instrument aberrations and forms a 
pattern of intensity variations or "speckle noise" on the 



final image. Part of the speckles cannot be calibrated 
because they evolve on various time scales (neither fast 
enough to smooth down to a halo nor stable enough to 
remove) and for this reason, these "quasi-static speckles" 
are one of the main limitations for high-contrast imaging. 

Several authors have discussed the challenge posed by 
the elimination of speckle noise in high-contrast multispec- 
tral images. It can be done by post-processing, after the 
best possible observations. Because images are highly spec- 
trally correlated, one can use the wavelength dependence 
of the speckles to subtract them. In the particular case of 
coronagraphic multispectral imaging, only some empirical 
methods have been developed to subtract the speckle field 
from the image in the focal plane. 

We propose an alternative approach based on a 
parametrized imaging model for the post-processing of mul- 
tispectral coronagraphic imaging corrected by an extreme 
AO system in the near-infrared domain. The aberrations 
and bright companions at small separations are estimated 
jointly in a Bayesian framework. In particular, it is possible 
to take advantage of prior information such as a knowledge 
on the aberration levels. This type of approach will be all 
the more efficient as the instruments improve with lower or 
more stable aberrations and more efficient coronagraphs. 

In section 2, we explain how previous methods used the 
information redundancy to suppress the speckles in high- 
contrast imaging. Then, we describe the advantages of a 
joint Bayesian estimation of the aberrations in the pupil 



1 



Ygouf et al.: Exoplanet detection and instrument aberration retrieval in multispectral coronagraphic imaging 



plane and of the planet map, based on a parametrized 
model of coronagraphic imaging. Section 3 presents the 
long-exposure coronagraphic imaging model that is used 
to simulate the images and to restore them. We also study 
the case of an approximate model. Section 4 describes the 
proposed Bayesian joint estimation method as well as the- 
oretical and numerical problems of an alternating restora- 
tion algorithm. In particular, we address the strong mini- 
mization difficulties associated to the aberration estimation 
and we propose some solutions. In Section 5, our method 
is validated by restoring images simulated with a perfect 
coronagraph. 

2. Post-processing speckle subtraction and 
multispectral imaging 

Several empirical post-processing methods have already 
been proposed to overcome the problem of detection lim- 
itation caused by the quasi-static speckles. Some of these 
methods used the wavelength dependence of the speckle 
pattern (Fig. [1} to estimate it and subtract it from the im- 
age, while preserving both the flux and spectrum of the 

pla net. 

I Racine et al.l (|1999() suggested to subtract two images 
at different wavelengths to eliminate both the point-spread 
function (PSF) and the speckle field in non-coronagraphic 
images. The main limitation of this simultaneous differen- 
tial imaging (SDI) method comes from the residuals caused 
by the evolution of the general PSF profile and of the 
speckle pattern with wavelength. These residuals can be 
reduced by increasing th e number of i mages used for the 
speckle field subtraction. iMarois etT al. ( 2000) showed with 
their double difference method that adding another image 
to the SDI theoretically improves the signal-to-noise ratio 
in the final image of the restored compani on. The case of 
multis pectral images has been tackled by ISparks fe Ford! 
(2002]), who described the so-called spectral deconvolution 
method in the framework of space-based observations for 
an instrument combining a coronagraph and an integral- 
field spectrometer (IFS). The method, subsequently im- 
pr oved and teste d on gr ound-based non-coronagraphic data 
by iThatte et al.l ()2007t ). is entirely based on a speckle in- 
tensity fit by low-order polynomials as a function of wave- 
length in the focal plane. More recently, Crepp et al. (2011) 
combined this method with the LOCI alg orithm, which is 
based on a linear combination of images (jLafreniere et al.l 
l2007fl . They tested this approach to restored on-sky images 
from the Project 1640 IFS on the Palomar telescope. These 
methods are applicable to any optical system and in partic- 
ular to those with coronagraphs. However, it is challenging 
to prevent the planet signals from being eliminated with 
the speckles because the planet presence is not explicitly 
modeled. 

In addition, some information on the measurement sys- 
tem can be very useful to dist i nguish a planet from the 
speckle field. iBurke fe Devanevl (|2010t ) combined classical 
empirical techniques of differential imaging with a multi- 
wavelength phase retrieval method to estimate the aberra- 
tion pattern in the pupil plane with a simple imaging model 
without a coronagraph. This multi-wav elength phase re - 
trieval is nicknamed wavelength diversity (|Gonsalvesll 982') . 
Information diversity is obtained here by different wave- 
lengths whereas it is obtained by introducing a known 
phase, e.g. defocus, in phase diversity. But in contrast to the 



phase diversity, the wavelen gth diversity does not re move 
the phase sign ambiguity. In lBurke fe Devanevl (|2010l ). the 
inversion algorithm is based on a maximum-likelihood es- 
timator, which measures the discrepancy between the data 
and an imaging model. The minimization of this estima- 
tor is all the more difficult as the number of unknowns to 
estimate is high. This problem is overcome by the sparse 
parametrization of the unknown phases 4>\ through the 
optical-path-errors (or aberrations) 5, assuming that the 
former are achromatic: </> (A) = 2ir5/\. This allows one to 
exploit jointly the images at all wavelengths to estimate 
the aberrations efficiently: the map of the unknown optical- 
path-errors 5 is common to all wavelengths. The number of 
unknowns is thus limited and the problem constrained. In 
the present case, Burke's wavelength diversity method does 
not apply readily, because it assumes non-coronagraphic 
imaging, whereas we consider the highly non-linear case of 
a coronagraphic imaging model. 

That is why we propose to take advantage of a com- 
bined use of wavelength diversity applied in a case of a 
coronagraphic imaging model, and a Bayesian inversion to 
jointly estimate the aberrations in the pupil plane and the 
planet map. The joint estimation aims at taking up the 
challenge of preserving the planets signal. An advantage of 
the Bayesian inversion is that it can potentially include an 
important regularization diversity to constrain the prob- 
lem, using for example prior information on the noise, the 
planet map (position, spectrum, etc.) or the aberrations. In 
the Bayesian framework, the criterion to be minimized is 
the sum of two terms: the data fidelity term, which mea- 
sures the distance between the data and the imaging model, 
and one or more penalty terms. An important difficulty is 
to define a realistic coronagraphic imaging model, that de- 
pends on parameters (e.g. aberrations) that can be either 
calibrated beforehand or estimated from the data. 



3. Parametric model for multispectral 
coronagraphic imaging 

To carry out the Bayesian inversion, we need a parametric 
direct model of coronagraphic imaging. This direct model 
will also be useful to simulate our test images. 

We used a non-linear analyt ical expression of the coron- 
agraphic image as proposed by ISauvage et al.l (|2010[) , with 
an explicit role of the optical aberrations before and after 
the coronagraph, and turbulence residuals. This model as- 
sumes that the coronagraph is "perfect" in the sense that 
the coherent energy is perfectly canceled out. The presence 
of upstream aberration however, will result in remaining 
intensity from the star in the image. The aberrations, or 
optical-path-errors, S, are assumed to be achromatic as an 
approximation. The most recent spectro-imagers take in- 
creasing care to avoid any source of chromatism, such as 
out-of-pupil aberrations, down to a level compatible with 
contrasts higher than 10 6 . The variable a = (a x , a y ) repre- 
sents the angular position in the focal plane in radians and 
the variable p = {p x , p y ) is the angular position in the pupil 
plane in radians -1 . Finally, \p = (Xp x , Xp y ) corresponds to 
a spatial position in the pupil plane in meters. 

We recall and discuss this model below. In particular, 
we estimate its simplified expression in the asymptotic case 
of very low phase, with its second-order Taylor expansion. 
This simplified expression helps to understand the explicit 
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950 nm I 1306 nm I 1647 nm 




Fig. 1. Evolution of the speckle field with the wavelength. Simulated images at 950, 1306, and 1647 nm for a 10 stellar 
flux over planet flux contrast. The dynamic range is adapted to the visualization. The speckle field moves with the wavelength but 
not with the planet position. 



way in which each type of aberration impacts the image. It 
also helps to identify some important ambiguities with dif- 
ferent sets of phases that can produce similar images, which 
will guide the subsequently selected approach for phase re- 
trieval. We also estimate the departure from this low-phase 
approximation when the phase grows and discuss the va- 
lidity of this approximation in a SPHERE-like case. 

3.1. Imaging model 

We assume that for an AO-corrected coronagraphic image 
at the wavelength A, the direct model is the following sum 
of three terms, separating the residual coronagraphic stel- 
lar halo, the circumstellar source (for which the impact of 
coronagraph is neglected), and noise ny. 

*a (a) = fx • h% (a) + [o x * hT] (a) + n x (a) , (1) 

where ix (a) is the data, i.e., the image to which we have ac- 
cess, /jj is the star flux at wavelength lambda and h^ c (a), 
the non-coronagraphic PSF, which can be estimated sepa- 
rately. Solving the inverse problem is finding the unknowns, 
namely the object ox (a) and the speckle field h\ (a), which 
we also call the "coronagraphic PSF" . 

3.2. Long-exposure coronagraphic PSF model 

A model description of h x (a) directly depends on the tur- 
bulence residuals and optical wave- front errors . After previ- 
ous w orks to model non-coronagrap hic PSFs (iPerrin et alJ 
200l and coronagr ap hie PSFs (ICavarroc et al.1 120061; 
Soummer et~aT1 l2007f ) . ISauvaee et alJ (|201Cj) proposed an 
analytical expression for the coronagraphic image with a 
distinction between upstream and downstream aberrations 
(cf. Equation (|16[) in Appendix B). The considered optical 
system is composed of a telescope, a perfect coronagraph, 
and a detector plane (cf. Figure @). Some residual tur- 
bulent aberrations 5 r (p, t) are introduced in the telescope 
pupil plane. 5 r (p, t) is assumed to be temporally zero-mean, 
stationary and ergodic. Because we only consider exposure 
times that are long with respect to turbulence timescales, 
these turbulent aberrations contribute only through their 
statistical spatial properties: power spectral density Sg r (a) 
or structure function D^. The static aberrations are sep- 
arated into two contributions: the aberrations upstream of 
the coronagraph S u (p), in the telescope pupil plane V u (p) 
and the aberrations downstream of the coronagraph Sd{p) 



in the Lyot Stop pupil plane Vd{p). The perfect corona- 
graph is defined as an optical device that subtracts a cen- 
tered Airy pattern of maximal energy from the electromag- 
netic field. Finally, the coronagraphic PSF depends on three 
parameters that define our system: the aberration maps 5 U , 
5d and, the residual phase structure function D^. 

3.2.1. Approximate long-exposure coronagraphic model in 
the low-phase regime 

Because the analytical expression for h\ is a highly non- 
linear function of the aberrations (jSauvage et al.|[2010l ). we 
derived and studied the relevance of an a pproximate ex- 
pression for this model (jYeouf et al.ll2010T ). Approximate 
coronagrap hic imaging models h ave been derived in sev- 
eral works. ICavarroc et al.1 (|2006H have developed a short- 
exposure expression and showed by simulations that the 
main limitation comes from the static aberrations and 
particularly the aberrations upstream of the coronagraph. 
Here, we consider a long-exposure imaging model and con- 
firm an alytically the dominan ce of the upstream aber- 
rations. ISoummer et al.1 (|2007l ) have developed a two- 
term expression with one static and one turbulent term. 
Nevertheless, these terms are not explicitly linked to the 
aberrations, which is what we are interested in. 

Assuming that all phases are small and that the spatial 
means of 4> u (p) and <frd(p) are equal to zero on the aperture, 
w e derive a second-orde r Taylor expansion of expression 24 
of ISauvaee et all (|2010h : 



[h c > 



, , 27T 

T 

/2tt 



V d (\p)*5 u (\p) 



Vd(Xp) *S Sr (a) 



P[S r (Xp,t)} 



V d (Xp) 



(2) 



where P d (Xp) and S u (\p) are the Fourier transforms of 
the downstream pupil and upstream aberrations and 
P [8 r (Xp, t)] denotes the piston of the aberration map 5 r (Xp, t) . 

(|P [S r (Xp, t)]\ 2 ) t ■ \pd ( Ap) j I is a corrective term that 

compensates for the fact that S r (Xp, t) is stationary and 
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Telescope pupil plane 
7^ : upstream pupil 



Lyot Stop 
pupil plane 
'Pd- downstream pupil 



Coronagraphic mask 



Detector focal plane 




<5^: static downstream aberrations 

S u : static upstream aberrations 

S r : residual turbulent aberrations with structure function Dg and power spectral density S$ 



°A * : circumstellar source 

■ : stellar coronagraphic halo 



Fig. 2. Optical scheme of a coronagraphic imager. The upstream and downstream static aberrations, and the adopted 
notations are denoted S u and Sd- Ai(a) denote focal plane complex amplitudes, whereas Vti(p) denotes pupil plane amplitudes. 



thus non-piston-free on the aperture at every instant. 3.2.2. Discussion 



\Vd (Ap) j is the Airy pattern formed by the pupil Vd (V) . 

This approximate expression brings physical insight into 
the long-exposure coronagraphic PSF model of Sauvage et 
al.: 

— The speckle pattern scales radially in A within the ap- 
proximate model and evolves as 1/A 2 in intensity in 
the data cube . It is consistent with the analysis of 
iSparks fc Fordl (|2002T ) . who performed fits of low-order 
polynomials as a function of the wavelength after rescal- 
ing radially. 

— The approximate expression can be separated into one 
static and one turbulent term. This is consistent with 
the analysis of ISoummer et ail (|2007f) with the advan- 
tage that these terms depend on the parameters of in- 
terest. The turbulent term is simply the turbulent aber- 
ration power spectral density, as seen at the resolution 
of the instrument, i.e., convolved by the output pupil 
Airy pattern. The static term is explicitly a function of 
the upstream aberrations. 

— The downstream aberrations do not appear in the static 
term. This confirms that the role of the aberrations up- 
stream and downstream of the coronagraph is very dif- 
ferent and that upstream aberrations are dominant in 
the final image. 

— Four equivalent upstream aberration sets, S u (p), 
fiu(-p), —fiu(p) and —S u (—p), which we call quasi- 
equivalent aberration maps in the following, lead to the 
same image (cf. Appendix A). This item is discussed in 
more detail in section R.4. 21 

— By using this approximate expression for h\ in the imag- 
ing model ([TJ), we can see that there is a degeneracy 
between the value of the star flux and the rms value of 
the aberration map, if there is no turbulence. Indeed, 
without turbulent aberrations, the approximate model 
multiplied by the star flux can be written as 



f*x-h c x (s u ) = r x 



(2nY 
A 2 



Pd*S u 



(3) 



Because of the complexity of the long-exposure corona- 
graphic PSF model of Sauvage et al., we first considered 
using the approximate model in our inversion algorithm to 
decrease the number of unknowns to estimate and simplify 
the criterion to minimize. But a study of this approximate 
model showed that the resulting image is too different from 
the one simulated with the Sauvage et al. expression: com- 
puting the root mean square of the difference between the 
two images leads to a sub stantial error of ty pically 30% in 
SPHERE-like conditions (|Ygouf et al.ll2010h . 

Speckle non-centrosymmetry. A substantial part of this er- 
ror arises because in the approximate model, the quasi- 
static speckles are purely centrosymmetric. Thus, if we 
eliminate the centrosymmetric part of the image by combin- 
ing it with a 180 degree rotated version of itself as follows: 



^antisym — 



*180 



(4) 



This is discussed in section 14.4.11 in greater depth. 



there are no residuals with this model, i.e. iantisym = 0. But 
this is not the case with a more physically realistic image 
such as the one simulated with the model of Sauvage et al.. 
The level of residuals after such a subtraction is determined 
by the quantity of upstream aberrations, as we can see in 
Figure |31 For example, with 30 nm rms of upstream aber- 
rations, the level of residuals is about six times lower than 
the rms value of the simulated image (cf. Figure 2]). Thus, 
for a significant quantity of upstream aberrations, using the 
model of Sauvage et al. rather than the approximate model 
for an inversion, should lead to fewer residuals on the final 
image. 



Speckle dilation. Another difference between the two mod- 
els tips the balance toward the model of Sauvage et al.. 
Indeed, in the approximate model the speckle dilation is 
powered by the 1/A 2 factor. If we subtract an image at 
950 nm from its 1650nm-dilation, there are no residuals 
with this model. The same operation with the model of 
Sauvage et al. (cf. Figure [5]) shows some residuals, at a 
level 2.5 lower than the rms value of the simulated im- 
age, which attests that it is not a pure speckle dilation. 
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To do this, we took the images at the minimum and max- 
imum wavelengths and rescaled the image at 950 nm with 
respect to the image at 1650 nm. Finally, we performed 
the following spectral differences between the two images: 
«diff 1650 = *i650nm - 7*950 nm, where 7 is a coefficient that 



minimiz es the squared diffe rence \i r _ 
given bv lCornia et all (|2010D : 



7*r. 



7 



Sp*950nm(/5)j 



1650 nm 



(P) 



J2p *1650nm(p) 



and is 



(5) 



A fine model of the speckle field must be able to account 
for deviations form centrosymmctry and for deviations from 
a radial scaling of the speckles proportional to the wave- 
length. That is why we adopted the model of Sauvage et al. 
rather than the approximate one in the inversion. 

3.2.3. Assumptions on the long-exposure coronagraphic PSF 
model 

The information we obtained from the approximate model 
study helped us define some key assumptions for the success 
of the speckle field estimation with the Sauvage et al. long- 
exposure coronagraphic PSF model. 

Because they have quite a different impact on the fi- 
nal image, it is important to distinguish the aberrations 
upstream and downstream of the coronagraph. The effect 
of the downstream aberrations is lower than that of the 
upstream aberrations, and furthermore, in predicted sys- 
tems such as SPHERE, they are expected to be much more 
stable and easier to calibrate than upstream aberrations. 
Additionally, because we consider long-exposure images, 
the residual turbulent aberrations will be averaged to form 
a smooth halo that is easily distinguishable from a planet. 
Furthermore, the statistical quantity D ( f >r , which charac- 
terizes this halo, will be measured t hrough the AO system 
wavefront sensor (|Veran et al.| [T997). Therefore, we here as- 
sumed that both the static downstream aberrations and the 
residual turbulent aberrations are calibrated and known. 
This decreases the number of unknowns because the only 
aberration map to estimate in order to access the corona- 
graphic PSF is the quasi-static upstream aberrations. We 
therefore denote the long-exposure coronagraphic PSF by 
h A {5 u ;5d,D,p r ) instead of h c x (S u , 8d, £ ) <^ r ) to underline the 
fact that 5d and D<f, r are assumed to be known. 

An advantage of our approach is that these assump- 
tions can evolve. The formalism will allow us to refine our 
method if we finally decide to estimate either the down- 
stream aberrations or the residual turbulent aberrations. 
Thus, we can slowly increase the complexity of the prob- 
lem in anticipation of using real data from SPHERE or 
from another instrument. 

4. Joint estimation of wavefront and object 
algorithm and minimization strategy 

This section introduces the criterion to be minimized (|4.1[) 
as well as the regularization elements that were used to 
constrain the problem for the present validations (|4.2I) . 
The minimization algorithm is then described, stressing 
the two stages that constitute its core (|4.3p . One of these 
stages presents some convergence difficulties. A minimiza- 
tion strategy is described in Subsection (|4.4I) . The chosen 
optimizer in described in Subsection 14.51 



4.1. Definition of the criterion to be minimized 

Following the Bayesian inverse problem approach, solving 
the inverse problem consists of finding the unknowns, firstly 
the object characteristics o (a, A) = {o\ (a)} A , secondly the 
parameters of the speckle field h x (8 u ; Sd, D^, r ) and /* (A) = 
{fx\\i which are the most likely given the data and our 
prior information about the unknowns. This can be reduced 
to minimizing the following criterion: 



j(oj*,s u ) = j2J2 



i 



ox*hZ c (S u ;6 d ,Dt r )\ 2 (a) 
R + R f , +Rs^ . 



\h - fx ■ h c x(S u ;8d,Dj, r ) 



(0) 



This criterion is the sum of two terms: the data fidelity 
term, which measures the distance between the data and 
the imaging model, and a non-exhaustive list of regular- 
ization terms on our unknowns R , R$. We consider 
that the noise is the sum of a photonic contribution and a 
detector contribution, that it is white and approximately 
Gaussian, which is a valid approximation for the flux as 
considered in this application. The noise variance is as- 
sumed to be known here and if it were not, it could b e 
estimated as a 2 n A = a^ h A + o\ et A (|Mugnier et al.|[2fJ04T) . 

where A = max(i,\,0) is the photon noise variance and 
o 2 det x is the detector noise variance previously calibrated, 
ft is assumed that the noise is not correlated from pixel to 
pixel or between images. 

The star flux at each wavelength can be analytically 
estimated from the criterion provided the regularization on 
flux is quadratic or absent. In the latter case, the maximum 
likelihood solution being given by Jpr = 0, we obtain: 



fx (o\,6 u ) 



y 




-ox* h n x c ) 




(a) 






(a) 



(7) 



Thanks to this analytical expression, the criterion to be 
minimized is that of Eq. ^ with / A replaced by / A , which 
will be denoted by J' (o,S u ) and depends explicitly on ox 
and S u only. 

Estimating both the object and the aberrations from 
a single image is a highly underdetermined problem. Any 
diffraction-sized feature in the halo can be interpreted ei- 
ther as a circumstellar point-source or as a part of the stel- 
lar speckle halo. This ambiguity can be decreased by using 
multispectral images but it is not sufficient. It is thus neces- 
sary to regularize the problem by adding more constraints. 
This is the role of the regularization terms. In this paper, 
we study the case of constraints on the object that is suffi- 
cient for the simulated images we used, but we should keep 
in mind that it is also possible to use constraints on the 
aberration map. 

4.2. Regularization terms and constraints 

We describe below the regularization terms that we used 
for the validation tests of this paper. Other regularizations 
could be chosen depending on the kind of images to be 
processed. 
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E 0.2 — 




upstream aberrations (n 



Fig. 3. Speckle non-centrosymmetry. Evolution of the level of residuals with respect to the rms value of the upstream quasi- 
static aberrations in an image simulated with the model of Sauvage et al. These residuals correspond to the antisymmetric part of 
the image. 




Fig. 4. Speckle non-centrosymmetry. In the same dynamic range: (left) simulated image, (center) 180 degree rotated version 
of the image and (right) residuals after combination of the image with the 180 degree rotated version of itself. The gain on the 
rms value after the subtraction is about 6. 




Fig. 5. Speckle dilation. In the same dynamic range: (left) image at 950 nm rescaled at 1650 nm and multiplied by the 7 
coefficient, image at 1650 nm (center) before and (right) after speckle subtraction. The gain on the rms value after the subtraction 
is 2.5. 



4.2.1. Regularization on the object R a 

A regularization on the object is fundamental to help com- 
pensating for the degeneracy that exists in the inversion 
between the aberrations and the object. By penalizing the 
energy in the object map, it favors the energy in the aber- 



ration map and prevents the speckles from being mistaken 
for a planet. 

The regularization term R Q includes the prior spatial 
and spectral information we have on the object. We chose 
here an L1-L2 white spatial regul arization, which assu mes 
independence between the pixels (jMeimon et all 12009ft be- 
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cause we are mainly looking for point sources. We used an 
L1-L2 regularization rather than a true LI regularization 
to keep a differentiable criterion, which simplifies the min- 
imization problem. The spectral prior is based on smooth- 
ness of the object spectrum. We currently assume that the 
object is white (constant spectrum) but as the final aim is 
to extract some spectra, for future va lidations we will use a 
L2 co rrelated spectral regularization (Thiebaut & Mugnier 
2006), which will involve the differences between the spec- 
trum at neighboring wavelength at each pixel and will en- 
force smoothness on the object spectrum. 

The regularization on the object is necessary to obtain 
a sparse object. Without the regularization many residuals 
remain on the estimated object. 

4.2.2. Positivity and support constraints on the object 

The object intensity map is a set of positive values, which 
is an important prior information. One should therefore en- 
force a positivity constraint on the object. This constraint 
can be implemented in various ways, such as criterion min- 
imization under the positivity constraint, reparameteriza- 
tion of the object, or explicit modification of the a priori 
probability dis t ributi on (e.g., addition of an entropic term). 
iMugnier et al.l (|2004[) have found that the best way to en- 
sure positivity, with respect both to speed and to not intro- 
ducing local minima, is to directly minimize the criterion 
under this constraint. We proceeded similarly. 

Because the star light is concentrated around the opti- 
cal axis, the flux is essentially estimated on this very bright 
region (cf. Equation [7]). But if there are too many residuals 
on the object, the flux estimation can be biased. Thus, im- 
posing, as is physically meaningful, that the object is null 
very close to the star, in a region of typically 3\/D radius, 
helps us estimate the star flux accurately. 

4.2.3. Regularization on the star flux Rf* 

To make the minimization more robust, it can be useful to 
constrain the flux estimation to physical values. Indeed, if 
there are no turbulent aberrations, there is an indetermina- 
tion between the star flux value and the phase rms value in 
the approximate model of Eq. (flU)) . Consequently, for aber- 
rations with very low rms values, the coronagraphic PSF h c x 
is close to zero and thus the analytical flux estimate given 
by Eq. can diverge. The presence of known turbulent 
aberrations naturally constrains the flux value, but to make 
the method more robust, we chose to prevent the flux from 
diverging. In practice, we regularized the flux estimation by 

(f*_f\2 

the following quadratic metric: Rf* = g 4 ■ This met- 

ric can be interpreted as a Gaussian prior low on the flux, 
but its role is not as essential for the criterion minimiza- 
tion as that of the L1-L2 regularization and the positivity 
constraint. This leads to the following expression for the 
analytic star flux: 

fx = o — • (8) 

E q (^a) 2 /<a + 1/^,a 

In practice, we chose a very high standard deviation <Jf t \ = 
100 x J2 a i\, to avoid biasing the flux. With this standard 
deviation, we can choose any mean flux, for example, /q = 



0. This is sufficient to avoid the division by zero in the flux 
computation and thus the flux divergence. 

4.3. Iterative algorithm 

The structure of the joint criterion of Eq. ([6]) prompted 
us to adopt an estimation of wavefront and object that 
alternates between estimation of the aberrations, assuming 
that the object is known (multispectral phase retrieval) and 
estimation of the object assuming that the aberrations are 
known (multispectral deconvolution) . 

Multispectral deconvolution. For given aberrations, we de- 
fine the following intermediate data where the (assumed 
known) stellar halo is subtracted, keeping then only the 
circumstellar object as seen in classical imaging: i' x ' = 
ix — f\ ■ h c x (5 u ; Sd, D^). By inserting these intermediate 
data into the the criterion of Eq. © , we obtain: 

J"(o, /*,*„) = ££ _ * J i'l - ox * hT(S u ;S d , Z^,.)| 2 (a) 

+ R + R f * +R 5 + ---, (9) 

which shows that the problem at hand is a non-myopic 
multispectral deconvolution of images i". The chosen regu- 
larization leads to a convex criterion (jMuenier et al.ll2004D 
and thus to a unique solution for a given set of aberrations 
and a given object regularization. 

Phase retrieval. If we replace the intermediate data i' x = 
ix — 0\ * h x lc (S u ;Sd,D^ r ) into the the criterion of Eq. ©, 
we obtain: 

J'(o, rA) = ££ - 5 1 -, J *' A - /a* • h c x (S u ; S d , D 4r )\ 2 (a) 

X Q Z(J n,X\ a > 

+ R + R f . +Rs + ---, (10) 

which shows that the problem at hand is essentially a phase 
retrieval problem. In this phase retrieval stage, the combi- 
nation of a high number of parameters to estimate (typi- 
cally 10 3 , see Section[S|) and of a highly non-convex criterion 
complicates the problem. To avoid local minima, several nu- 
merical solutions resulting from a fine understanding of the 
imaging process are necessary and are described below. 

4.4. Phase retrieval: dealing with local minima 

4.4.1. Choice of an appropriate starting point: very small 
random phase 

To keep the computation time reasonable, we used the local 
descent algorithm described in Subsection 14.51 to minimize 
the criterion. Because the latter is highly non-convex, the 
chosen starting point so that we are fully in the conditions 
where the Taylor expansion developed in 13.2.11 is valid and 
where the criterion is less non-convex. It allows the algo- 
rithm to avoid many wrong directions, and thus many local 
minima. As the algorithm converges, the upstream aberra- 
tion rms value increases toward its true value and a gradual 
non-linearity of the model is gradually introduced. 

We tested the phase retrieval capability of our algorithm 
with respect to the chosen starting point, assuming that 
there is no object to estimate. We give the rms value of 
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the estimated upstream aberration map and the rms value 
of the difference between the simulated and the estimated 
maps estimated as follows: 



rms 



diff 



simulated 



^estimatcd^ * 



1/2 



EM 



simulated^ 



1/2 



x 100. (11) 



The inversion is performed with one spectral channel and 
without turbulence. Figure H2 compares some estimated up- 
stream aberration maps (a, b, c) to the simulated one 
("true"): 

(a) Using a random aberration map with the same rms 
value as the true aberrations (30 nm at 950 nm) as a 
starting point does not help in finding the global mini- 
mum. Indeed, the algorithm converges very quickly to- 
ward a local minimum and the estimated aberration 
map (rms a = 307 nm at 950 nm) is completely differ- 
ent from the simulated one (rms dlff ' a = 10 3 %). 

(b) Using a zero-aberration map as a starting point does 
not work either. This is probably because the approx- 
imate model is an even function. For this particular 
starting point, the gradient is null, which leads to some 
convergence difficulties. The rms of the difference be- 
tween the two maps is about 1.4 x 10 4 %. The estimated 
aberration pattern (rms b = 4069 nm at 950 nm) seems 
to show that the algorithm does not explore the high 
frequencies. 

(c) The solution we propose is to use as a starting point 
for the minimization a non-null random aberration map 
with a low rms value compared to those of the "true" 
simulated aberration map. In practice, we chose an rms 
value about 10 8 times lower than the "true" value. This 
leads to a correct estimation of the aberration map 
(rms c = 30.2 nm at 950 nm) with an rms of the dif- 
ference between the two maps of about 0.6%. 

If we plot the same results for images simulated with tur- 
bulence, the conclusions are not the same. The convergence 
to an aberration map that resembles the true one and has 
similar rms value is easier: 

(a') By using a random aberration map with the same rms 
value as the true aberrations (30 nm at 950 nm) as a 
starting point, the rms value of the estimated aberration 
map (30.7 nm) is close to the true value. After a careful 
inspection, the estimated aberration map turned out to 
be similar to the opposite of the simulated one. This is 
discussed in more detail in 14.4.21 

(b') Using a zero-aberration map as a starting point leads 
to a good pattern and a good rms value of aberration 
map (30.1 nm). 

(c') Using a non-null random aberration map with a low 
rms value as a starting point also leads to a good pattern 
and a good rms value of the aberration map (30.1 nm). 

For the (b') and (c') cases, the difference between the two 
maps is bigger than before (17%) but the estimated aber- 
ration map is a sufficiently good starting point for the al- 
ternating minimization. 

The choice of an appropriate starting point seems not to 
be as essential with turbulent aberrations as it was with- 
out turbulent aberrations. Indeed, the presence of turbu- 
lent aberrations raises the ambiguity that exists between 



the value of the star flux value and the rms value of the 
upstream aberration map. Nevertheless, even if we assume 
that there are turbulent aberrations in the following, we 
chose to use a random aberration map with a low rms value 
as a starting point of the phase retrieval because it allows 
us to avoid some local minima by linearizing the highly 
non-linear model used in the inversion. 



4.4.2. Avoiding some local minima by testing 
quasi-equivalent starting points 

In the approximate model, four different aberration maps 
can give the same image (cf. Equation (|15p in Appendix 
A). This means that, from a given starting point, the min- 
imization algorithm can take four different but equivalent 
directions from the approximate model point of view. But 
from the point of view of the model used in the inversion, 
this is not the case because it depends on downstream aber- 
rations, which break the symmetry. That is why we call 
them "quasi-equivalent" aberrations maps (cf. 13.2.1]) . 

The idea is then to explore the several regions offered 
by the four different quasi-equivalent aberrations maps to 
determine which of these solutions gives the smallest crite- 
rion. To do this, we performed an initialization step where 
the very small random phase is taken as a starting point. A 
first phase retrieval stage was performed with this starting 
point, leading to a first estimated aberration map denoted 
by 5 u init,1 (p). Then, the three other quasi-equivalent aber- 
ration maps S u zmt,1 (— p), —S u tn1 ' t '(p) and — (5 M mlt,1 (— p) 
were taken as starting points for three other phase retrieval 
stages. This led to three more estimated aberration maps 
denoted by 5 u mlt ' 2 , 6 u mtt ' 3 and 5 u mtt < 4 . 

Figure (J7J shows the four estimated aberration maps at 
the end of the initialization step. These estimated aberra- 
tion maps are compared to the simulated one (Fig. ([BJd)). 
The final aberration map chosen as a starting point for the 
alternating algorithm is the one that gives the minimum 
value for criterion J of Equation H2 

(6 u ) lmt = arg min { J [$™ M ] , J ft""' 3 ] , J [C**' 3 ] , J [C UA ] } • 

(12) 

In practice, quite often and in particular in this simulation, 
it turns out that the chosen set of aberrations is also the 
one whose rms value (rms b = 30.1 nm at 950 nm) is closest 
to the 'true" phase (Fig. ([H "true") and Fig. (JBIfl "true")). 



4.4.3. Avoiding some local minima in the multispectral 
inversions by taking the previously estimated 
aberration map as starting point 

We used a local descent algorithm to minimize the criterion. 
For this, the gradients are computed for each explored di- 
rection and the computation is all the longer as there are 
spectral channels. That is why it is useful to perform several 
inversions by gradually increasing the number of spectral 
channels. This also helps us avoid some local minima. We 
begin by an inversion with one spectral channel. Then, we 
add some more spectral channels for new inversions and 
each time, we take the previous estimated aberration map 
as a starting point. 

Because an inversion with only one spectral channel 
sometimes leads to a local minimum, it is useful to also test 
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Fig. 6. Choice of an appropriate starting point. Estimated upstream aberration maps with one spectral channel for three 
different starting points, [i] Without turbulent aberrations in the simulated images, [ii] With turbulent aberrations in the simulated 
images. From left to right, with a dynamic range adapted to the visualization: "true" simulated aberration map, (a) and (a') 
estimated aberrations with a random aberration map (rms value of the simulated aberrations) as starting point, (b) and (b') 
estimated aberrations with a zero aberration map as starting point, and (c) and (c') estimated aberrations with a random aberration 
map (rms value 10 8 times lower than the "true" one) as starting point. The estimation is performed with a regularization on the 
star flux. 




rms 1 — 31.2nm rms 2 — 30.1 nm rms 3 — 31.7nm rms 4 — 31.6nm 

rms diff ' 1 = 160% rms diff ' 2 = 17% rms diff ' 3 = 140% rms dift ' 4 = 160% 



Fig. 7. Estimated upstream aberrations for the four quasi-equivalent aberration maps as starting points. From left 
to right, with the same dynamic range: , <5™*' 2 , i^ 7 " 4,3 , J™ 4,4 . The image simulation is performed with one spectral channel 

in the presence of turbulent aberrations. 



the four quasi-equivalent starting points with two spectral 
channels (cf. Section |4.4.2[) . 



4.5. VMLM optimizer 

To minimize the criterion, we chose the variable met- 
ric with limited memo ry and bounds (VMLM-B) 
method (|Thiebautl |2002j) Updated from the BFGS 
variable metric method (jPress et al.1 l2007t ). it is us- 
able for a problem of large dimensionality. Moreover, 
it offers the possibility to constrain these parameters. 
This makes this method a good tool for many inversion 
probl e ms in high an g ular resolution (jMeimon et al.l 
120091: iGratadour et all I2005Q. It is available from 
http : //www-obs .univ-lyonl . f r/labo/perso/ eric . thie 



4.6. Summary of the developed algorithm 

Figure © summarizes the different steps of the developed 
algorithm. The choice of a very small random phase as a 
starting point is essential because it avoids falling into some 
local minima (Section 14.4. An initialization phase is per- 
formed. It consists in running the algorithm for the four 
quasi-equivalent solutions (Section l4.4.2p . The solution that 
leads to the lowest criterion value is selected. Then, the 
minimization core is performed, alternating between the 
aberration estimation, assuming that the object is known 
(multispectral phase retrieval), and the object estimation, 
assuming that the aberrations are known {non-myopic mul- 
tispectral deconvolution) . Regularization terms and con- 
straints prevent the algorithm from falling into other local 
fflgfafirpfo i|[fc[3ig:tk)Jh te2.,3p . Iterations are performed (Section 
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14. 3p until the stopping rule of the algorithm is ver ified. The 
chosen optimizer is the VMLM-B (|Thiebaudl2002h (Section 
153]) . 



5. Validation of the inversion method by 
simulations 

In this section, we validate the exoplanet detection capa- 
bilities of our inversion method. After giving the numerical 
simulation conditions, we investigate the estimation quality 
of the aberrations and the object as a function of the num- 
ber of images at different wavelengths used. We also study 
the algorithm robustness with respect to the simulated im- 
ages and with respect to the starting point we use. Finally, 
we study the effect of the bandwidth on the quality of the 
object estimation. 

5.1. Simulation hypothesis 

From a data cube of six images simulated with the im- 
age fo rmation model of equation (fTJ) and the ISauvaee et al.l 
(2010) analytical expression for coronagraphic imaging (cf. 
Equation (|16p in Appendix B), we jointly estimated the 
speckle field and the object map. We chose pixel indica- 
tor functions as the basis for the phase rather than, e.g., 
a truncated basis of Zernike polynomials, to model and re- 
construct phases with a high spatial frequency content. The 
hypotheses are typical of a SPHERE-like instrument: up- 
stream 5 U and downstream Sd aberrations simulated with 
standard deviation of 30 nm (cf. their power spectral densi- 
ties in Figure (TlU)) ). star-planet angular separations of 0.2 
and 0.4 arcsec, contrasts, i.e. ratio of star flux over planet 
flux of 10 5 , 10 6 and 10 7 , a [950 nm; 1647 nm] spectral band- 
width and an integrated flux of 4 x 10 11 on the data cube 
in presence of photon noise and a transmission (throughput 
and quantum efficiency) of 10%, corresponding to the obser- 
vation of a 6-magnitude star for 25 minutes with the VLT. 
We used 128 x 128 pixels to simulate our images, Shannon- 
sampled at 950 nm. This results in a number of unknowns 
to estimate for the aberration map of about 3 x 10 3 . If wc 
add the unknowns to estimate for the object map, which is 
16 x 10 3 , the total number of unknowns is about 2 x 10 . 
Figure © shows the simulated objet map (|TJa]) and the 
associated image in the focal plane (|9b|) . For an easier vi- 
sualization, we represent the images in the focal plane and 
not the object map in the following. Figure [TJ] shows the 
simulated aberration map (|9"c)) and the associated image of 
the speckle field in the focal plane (|9dl) . 

5.2. Algorithm robustness and performance studies 

We jointly estimated the upstream quasi-static aberration 
map and the object map with multispectral data. To study 
the robustness of the method we have developed, we ran 
several simulations in a Monte Carlo-like manner. Both dif- 
ferent simulated images ([5.2. ip and different starting point 
(15.2.2j) were used for the inversion. The results of the inver- 
sions with two and six spectral channels were compared to 
study the effect of the redundancy of information offered 
by the multispectral imaging. 




a Object map b Image of the object map 

in the focal plane 




c Aberration map d Image of the speckle field 
in the focal plane 



Fig. 9. Simulated images at A = 950 nm. (a) Simulated ob- 
ject map and (b) associated image in the focal plane. The follow- 
ing planets are simulated: one with a star-over-planet contrast of 
10 at a separation of 0.2 arcsec (planet 1), two with star-over- 
planet contrast of 10 6 at separations of 0.2 (planet 2) and 0.4 
(planet 3), respectively, and one with a star-over-planet contrast 
of 10 7 at a separation of 0.4 arcsec (planet 4). The image in the 
focal plane is obtained by convolving the object map o\ by the 
non-coronagraphic psf h" c . (c) Simulated aberrations and (d) 
associated image of the speckle field in the image focal plane. 
The image is given by the coronagraphic PSF h^- 



5.2.1. Different simulated phases 

We applied our method to ten images simulated with dif- 
ferent random upstream aberration maps to assess the al- 
gorithm robustness. Figure QT] shows the estimated images 
of the object from these images for a two-spectral channel 
inversion (a) and a six-spectral channel inversion (b). These 
results bring several conclusions, both on the robustness of 
the method and the multispectral redundancy. 

Robustness of the method. In one out of ten cases, some 
residuals from the speckle field remain on the object. The 
level of this residual prevents one from detecting any planet. 
In the other nine cases, the level of the residuals is negligi- 
ble, which allows one to detect at least one planet. These 
observations, independent of the number of spectral chan- 
nels used for the inversion, shows that the method is rela- 
tively robust, given the minimization difficulties we met. 

Multispectral redundancy. We can see some significant dif- 
ferences between the two-spectral channel inversion and the 
six-spectral channel inversion: 
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Fig. 8. Block diagram of the algorithm used for the joint estimation of the object map and the upstream aberra- 
tions. 



— The planet with a contrast of 10 5 is detected in eight 
out of ten cases with two spectral channels and in nine 
out of ten cases with six spectral channels. 

— Only one planet with a contrast of 10 6 is detected with 
two spectral channels in three out often cases. The same 
planet is detected in nine out of ten cases with six spec- 
tral channels. 

— The other planet with a contrast of 10 6 is detected in 
three out of ten cases with six spectral channels. 

— Furthermore, the estimation of the flux of the planets is 
more accurate with six spectral channels. The inversion 
shown in Figure [T2] is representative of a large body 
of simulated tests and is performed with two and six 
spectral channels. The simulated and estimated image 
of the object maps for these two different inversions are 
represented as the errors on the estimated planet flux 
when the planet is effectively detected. 



These results show that multispectral redundancy helps us 
detect the planets. However, we note that even if the conver- 
gence problems were solved with our minimization strategy 
in most cases, there are still challenges to be overcome to 
arrive at a perfectly robust algorithm. 



5.2.2. Different starting points 

For each of the ten previous simulated phases, we took ten 
different random aberration maps with the same level of 
rms value as starting points to demonstrate than any ran- 
dom aberration map, provided it is small, allows us to find 
a good solution. Figure [T3] shows the estimated images of 
the object for the ten different random aberration maps, 
used as starting points, for one phase. 

The starting point does not have a strong impact on the 
planet detection. In all cases, the three brightest planets 
are detected. In only one out of ten false planet is 

detected. Despite all the attention given to the initialization 
of the algorithm, the choice of starting point can have, to 
a limited extent, an effect on detection. 

5.3. Bandwidth effect 

We study here the bandwidth effect on the planet 
image estimation quality. We selected the following 
spectral bandwidths: [950 nm; 1050 nm], [950nm;1150nm], 
[950 nm; 1350 nm], and [950 nm; 1650 nm]. The two last 
bandwidths correspond to operational modes of the 
SPHERE- IFS instrument, whereas the other bandwidths 
are produced for the purpose of our study's completeness. 
Figure (|14j) compares the estimated planet image maps for 
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Power Spectral Density of Upsiream Aberrations 



Power Spectral Density of Downstream Aberrations 




a Power spectral density of the upstream aber- b Power spectral density of the downstream 
rations. aberrations. 

Fig. 10. Power spectral densities of the simulated aberrations. The upstream and downstream aberrations are randomly 
generated according to a f~ 2 spectrum and with an rms value of 30 nm. The downstream aberrations are not corrected by the 
adaptive optics. The upstream aberrations are corrected by the adaptive optics up to 20 cycles/pupil. The residuals are due to 
the non-common path aberrations. Below 4 cycles/pupil, these non-common path aberrations are corrected but there are some 
residuals from rotative optics. 




a two-spectral channel inversion 



'■•■t^t 



b six-spectral channel inversion 

Fig. 11. Robustness study on the simulated phases. With the same dynamic range, at 950nm: estimated planet images 
[o\ * h^ c ] (x, y) from different images, simulated with ten different randomly generated upstream aberration maps. 



these different inversions. If the planet with a contrast of 
10 5 is detected each time, only the two broader bandwidths 
(700 and 400 nm) allow one to detect the two planets with 
a contrast of 10 6 . With a 200 and a 100 nm bandwidth, the 
planet with a contrast of 10 6 which is at a separation of 



0.4arcsec is detected, whereas the one that is at a separa- 
tion of 0.2arcsec is not. 

The detection performance increases with the band- 
width because the incorporated spectral information helps 
the algorithm to distinguish between speckles and planets. 
Indeed, from one wavelength to another, the amplitude of a 
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flux error planet 1 = 19% 
flux error planet 3 = 82% 



flux error planet 1 =1% 
flux error planet 2 = 73% 
flux error planet 3 = 47% 




Fig. 12. Multispectral redundancy. With the same dynamic range, at 950nm, from left to right: simulated planet image 
[o\ * h^ c ] (x,y) and estimated ones with a two- and six-spectral channel inversion. The errors on the planet flux estimations are 
computed when the planet is detected. 



Fig. 13. Robustness study on the aberrations maps used as starting points. With the same dynamic range, at 950nm: 
estimated planet images [o\ * /i^ c ] (x, y) from one image with different random aberration maps used as starting points. 



speckle movement is proportional to the wavelength differ- 
ence and to the radial position of the speckle. For a given 
spectral bandwidth, it is therefore easier to detect a planet 
that is far away from the star than a planet than is close 
to the star. 



6. Conclusion 

We have proposed an original method that jointly estimates 
the object (multispectral deconvolution) and the aberra- 
tions (multispectral phase retrieval) for the new genera- 
tion of planet finders. For the first time, a fine parametric 
model of coronagraphic imaging, describing the instrument 
response, is used for the inversion of simulated multispec- 
tral images in a solid statistical framework. Even though 
the model remains a simplification of reality, in particular 
when assuming achromatic wavefront errors, it goes much 
further than only assuming that the spatial speckle pattern 
essentially scales with wavelength. We have shown that the 
second-order approximation of the imaging model has the 
same behavior as the often-used model for the problem at 
hand: the speckle pattern is centrosymmetric, it scales with 
the wavelength, and it is not dependent on the downstream 
aberrations. The departure from this case quickly becomes 
very significant as the phase grows. It is clear that the abil- 
ity of a finite phase to produce non-symmetric speckles 
induces a strong ambiguity between the estimates of the 



aberrations and the object, even if the latter is a point-like 
companion. 

To set up our method, we developed an iterative al- 
gorithm. With only one spectral channel, this joint esti- 
mation is an underdetermined problem, as we emphasized 
before. This underdetermination results mathematically in 
a degeneracy of the global minimum. A multispectral inver- 
sion raises this underdetermination but it is still possible to 
fall into local minima. Because of the high non-linearity of 
the coronagraphic imaging analytical model and the num- 
ber of unknowns to estimate (about 10 3 in our case), the 
phase retrieval, even if it is multispectral, remains a difficult 
problem. We set-up a minimization approach that remains 
quite fast (without systematically exploring the whole pa- 
rameter space in a search for global minimum) but that is 
still relatively robust: extensive tests showed the success in 
converging to a good solution in 90% of the cases in a sys- 
tematic manner, and in the other cases, the failure appears 
obviously in the results with no risk of confusion and can 
motivate tests with alternative criterion minimization ap- 
proaches. We obtained these convergence capabilities of the 
algorithm by bringing original solutions to the minimiza- 
tion difficulties of the phase retrieval, inspired by studying 
the imaging model. One element of the solution is to use 
a small random aberration map as starting point. Another 
element is to explore the several directions offered by the 
quasi-equivalent aberration maps. 
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Fig. 14. Bandwidth effect. With the same dynamic range, at 950nm, from left to right: estimated planet images [o\ * h^ c ] (x, y) 
with a 700, 400, 200, and 100 nm-bandwidth inversion. 



A wide variety of prior information, either about the 
system (aberrations, flux, noise) or about the object of in- 
terest, can be used to constrain the problem. The choice of 
a Bayesian approach allows this flexibility. In particular, a 
regularization that minimizes the energy on the object map 
helped us in separating the aberrations from the object and 
in decreasing the speckle noise in the reconstructed object 
map. 

The restoration of images simulated with a perfect coro- 
nagraph is very encouraging for the extraction of planetary 
signals at levels that begin to be astrophysically interest- 
ing. We demonstrated the efficiency of the method even 
with only two spectral channels, by achieving a contrast of 
10 5 at 0.2arcsec. Multispectral redundancy improves the 
detection, which allowed us to achieve a contrast of 10 6 at 
0.2arcsec with six spectral channels. 

We therefore believe that this approach will be quite 
powerful when we are faced with experimental data. This 
deserves to be studied, as well as how the performance 
will evolve in the cases of images simulated with a non- 
perfect coronagraph, real images from the SPHERE instru- 
ment in the lab, or real images from an instrument on-sky. 
Eventually, this method could be used to improve the per- 
formance of the existing multispectral imaging instruments, 
providing better astrophysical exploitations. Now that we 
demonstrated that we can manage the difficulties linked to 
the criterion minimization, we can now focus on its applica- 
bility, and on adding more prior information, starting with 
the full set of information that can be obtained from the 
instrument calibrations. 

The lessons learned by applying the method could also 
facilitate the approach for the design of future instru- 
ments such as EPICS for the European Extremely Large 
Telescope ()Kasper et al.l [20081) . and the definition of their 
calibration procedures. 



without the variables for better readability: 

, 2 



Vd*(j>u 



+ 



v d 



*S<f, r (a)-{\P\ 
'■)■ 



Vd 



(13) 



We consider here the static term 



Pd*<Pu 



and re- write 



it in the form of a correlation. For this, we consider two 
functions f — Vd and g — 5 U of the two variables p x and 
p y and denote f(p) = f(-p) = f(-p x , -p y ). 

By using the definition of the correlation T f g (p) and the 
convolution Cf g (p) of the two functions f(p) and g(p), 



r /9 (p) = Hp) ® g{p) = J r(p')g(p' + p)d P , 

c fg(p) = f(p)*9(p) = J f(p')g{p - P'W, 
and the properties 

(/)* = J* andf*g = f®g*, 



we obtain 



f*9 



fg ■ (fg)* = fg ■ (fg)* = fg* (fg)* = fg®fg = T fg . 



This yields 



= r 



(14) 



The properties of the autocorrelation 



r 



/(p) 



r 



and 



Appendix A: Indetermination on the estimated 
aberrations from an image simulated with our ap- 
proximate model 

We show here that four sets of upstream aberrations give 
the same image for our approximate model. We re-write 
the expression below as a function of <j) u = (2tt/X) x 5 u and 



lead to 



and 



r -/(p) - r /(p)' 



T (V d -S u )(p) - r (7V<5„)(-p) 
T (VH-Su))(p) = T (Td-(Su))(-p) 



T (V d -(-S u ))(p) - r (7V<M(p)' 
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This yields 



CP d -S u )(p) 



(V d -5 u ){-p) 



(■PH-s u ))(p) 



r 



(15) 



which means that the upstream aberration sets S u (p), 
$u(—p), —$u(p) and ~&u(—p) are equivalent with respect 
to the approximate model, because they give the same im- 
age. This is true even in the presence of the turbulent term 

2^ 



Pd 



*V (a)-(\P[4> r }\ 



Appendix B: Indetermination on the estimated 
aberrations from an image simulated with the 
Sauvage et al. model 

In classical imaging, i.e. "non-coronagraphic imaging" , the 
sign of the even part of the phase can not be dedu ced 
from only one image in the focal plane ()Bland l2002h . In 
other words, if we denote 4> e and 4> , the even and the 
odd parts of the phase, the two phases (f> = <j) e + <p and 
4>' = —<j) e + <f) give the same image. In this appendix, we 
show that this is also the case for the Sauvage et al. expres- 
sion (Sauvage et al.ll20Tol) in coronagraphic imaging, if we 
assume that the sign of the even part changes for all phase 
errors. 

The expression of Sauvage et al. is 
h\ = (A n A* n ) + (\( Vo )\ 2 )A d A* d -2Vl{(r la A; i )A d }, (16) 

with A n \a) = TF _1 [V d (p) ^""W] , A d (a) = 

TF" 1 [V d (p) e^ d ^] , ct>i{p) = 2ir^ and <t> tot (p) = 
4> r (p)+4> u (p)+<j) d (p). TF [.] denotes the Fourier Transform. 

^|(?7o)| 2 ^ represents the mean Strehl ratio during observa- 
tion, such as 



T) (t) = (* (p)|n(p)) 
1 



%(p)V u (p)d 2 p 

Pl(p)e-^d 2 Pl 



S 
1 

TP 



with 4>(p,t) = (j> r {p,t) +<p u (p). 

The first term of Sauvage et al.'s expression (A„A*) 
is the classi cal case of non-coronagraphic PSF, which is 
well-known (jBland 120021) . The term A n A* n stays identical 
whatever the sign of the even part of the phase. 

The second term of Sauvage et al.'s expression is the 

product of two factors: ^|(?7o)| 2 ^ and A d A* d . The latter stays 

identical whatever the sign of the even part of the phase. 
We take the following phase <f>' = — <fi e + <\> Q and calculate 



the corresponding r]' (t), assuming that p" = —p: 
Vo(t) 



1 f[„2, \ 


'IV (P'l^p 


1 ft 9/ s 




i r r „ 
° j j p 


{<i>e(-P,t)+<t>o(-P,t)]d 2 p 




e j[0 c (p",t)+0 o (p",t)] d 2 






[»(«)]■■ 





K^o)! J = ( \(Vo ■ %(*)*) I ) is then independent of the 

sign of the even part of the phase. Thus, the product 

^|(?7o)| 2 ^) A d A* d is also independent of the sign of the even 

part of the phase. 

We study now the third term 25i {(r]oA^) A d }. Assuming 
that (j) d = (<j> d ) e + (<j) d ) and <fi' d = -((j> d ) e + {<j> d ) - 

A' d {a) = TF" 1 [V d {p)e*+'V>' 

V d (p) e^H&OeCpHOMoO)]] e~ 2 " {pa) A 2 p 
V d (p") e -J'[C^)e(p")+(^)<>(p")]l e 2m ( p " a \l 2 p" 



= [Ad]*- 

In the same way as the previous demonstration, we can 
show that A* = A n . Under the effect of the transformation 
4> — >• 4>' , the different terms become 




In other words, (rjoA^ ) A d -4 [{rj Q A* n ) A d ]* When we take 
the conjugate of a complex number, only the sign of the 
imaginary part changes. Because we take the real part of 
this expression, changing the sign of the even part of the 
phase does not change the term. 

To conclude, like in classical imaging, changing the sign 
of the even part of the phase does not change the image 
in the focal plane. This means that two sets of aberrations 
give the same image. But if we assume like in this com- 
munication that the downstream aberrations are fixed and 
known, this removes the degeneracy. 
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